Resonant generation of coherent phonons in a superconductor by ultrafast optical pump pulses 
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We study the generation of coherent phonons in a superconductor by ultrafast optical pump pulses. The 
nonequilibrium dynamics of the coupled Bogoliubov quasiparticle-phonon system after excitation with the 
pump pulse is analyzed by means of the density-matrix formalism with the phonons treated at a full quan- 
tum kinetic level. For ultrashort excitation pulses, the superconductor exhibits a nonadiabatic behavior in which 
the superconducting order parameter oscillates. We find that in this nonadiabatic regime the generation of co- 
herent phonons is resonantly enhanced when the frequency of the order-parameter oscillation is tuned to the 
phonon energy, a condition that can be achieved in experiments by varying the integrated pump pulse intensity. 
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I. INTRODUCTION 

The generation of coherent phonons by ultrashort optical 
pulses with duration much shorter than the phonon vibra- 
tional period has been extensively studied in various materi- 
als, such as bulk semiconductors^— semiconductor quantum 
wells^ii and superlatticesr^i£ as well as high-temperature 
superconductors^— For semiconducting systems, several 
distinct coherent phonon generation mechanisms have been 
discussed^^ For example, in the displacive mechanism, the 
optical pulse creates a finite photocarrier distribution almost 
instantaneously on the time scale of the phonon subsystem^ 
This results in an abrupt change of the equilibrium positions 
of the lattice ions, and hence gives rise to coherent oscilla- 
tions of the atoms around the new potential minima. In the 
impulsive mechanism, an effective direct coupling of the laser 
field to the lattice ions is assumedji^ 2 . leading to a brief and 
intense force acting on the atoms. The detailed dynamics of 
the electronic subsystem on time scales longer than the optical 
pulse is, in general, irrelevant for the description of coherent 
phonon creation in semiconductors. An exception to this rule 
occurs when the electronic subsystem oscillates with a period 
on the time scale of the phonon vibrations, in which case the 
coherent phonon generation is resonantly enhanced . This has 
been observed both in semiconductor quantum and 
in superlatticesji^i^ 

In this paper, we investigate the generation of coherent 
phonons in a nonequilibrium superconductor. Specifically, 
we study the optical excitation of Bogoliubov quasiparticle 
states above the superconducting ground state on time scales 
shorter than the phonon vibrational period r p h. We find that, 
akin to the displacive mechanism in semiconductors, a sud- 
den change in the Bogoliubov quasiparticle distribution func- 
tions generates coherent phonon oscillations. This mecha- 
nism of phonon creation is relevant to pump-probe experi- 
ments on superconductors with a pump photon energy of the 
same order but slightly larger than twice the superconduct- 
ing gap amplitude |A| and a laser pulse duration r p shorter 
than both the phonon period r p h and the dynamical time scale 
of the superconducting order parameter ta ~ h/(2\A\). It 



has recently been shown that, whenever the pump pulse du- 
ration t p is much shorter than ta, oscillations are created 
in the quasiparticle occupations with frequency of the order 
of ~ 27t/ta ~ 2|A|/7i<2 3 .— For r A < r ph , these oscilla- 
tions average out on the time scale of the phonons and are 
therefore unimportant for the creation of coherent phonons. 
When ta is close to r p h, on the other hand, the generation of 
coherent phonons is resonantly enhanced. Remarkably, pro- 
vided that o; p h < 2|A|/ft, the Bogoliubov quasiparticle oscil- 
lations can be brought into exact resonance with the phonon 
frequency u p h by adjusting the integrated pump pulse inten- 
sity (see Figs. Q] and |2). 



In the following, we theoretically investigate this resonant 
coherent phonon generation mechanism by employing a mi- 
croscopic model of an s-wave superconductor coupled to an 
optical phonon mode with frequency u; p h. We study the pulse- 
induced dynamics of this model system at times shorter than 
the quasiparticle energy relaxation time r e , a regime which 
can be fully described within mean-field BCS theory Dif- 
ferent orderings of the involved time scales are studied with 
a particular emphasis on the case where both the phonon and 
the quasiparticle subsystems evolve in a nonadiabatic fashion, 
i.e. where r p <C T p h , t a <C r e . In this nonadiabatic regime, 
traditional approaches for computing nonequilibrium dynam- 
ics, such as the time-dependent Ginzburg-Landau theory or 
the Boltzmann kinetic equation, are not applicable, since the 
full dynamics of both the normal and the anomalous quasi- 
particle densities, as well as that of the coherent-phonon am- 
plitudes needs to be accounted for. Therefore, we resort to 
the density-matrix formalism^ 4 , to numerically compute the 
coherent response of the model system after excitation by a 
short pump pulse. Based on this approach, we analyze in de- 
tail the generation of coherent phonons and calculate lattice 
displacements both for resonant and off-resonant conditions. 
The analysis presented in this paper is complementary to the 
one of Ref. |35[ which employs Boltzmann-type kinetic equa- 
tions to study the adiabatic dynamics of Bogoliubov quasipar- 
ticles coupled to incoherent phonons. 
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II. MICROSCOPIC MODEL 

The microscopic model we consider is a single-band BCS 
s-wave superconductor coupled to an external electromag- 
netic field and to a single branch of optical phonons H = 
H sc + H em + H p h + -ffe-ph- Within mean-field theory, the BCS 
superconductor is given by the following Hamiltonian 

h sc = £ kcL c k, CT + E [ Ac V-ki + A * c - k i c kt 

k,<r kew 

(la) 

where Ck.er represents the electron annihilation operator with 
spin a and momentum k, £k = fi, 2 k 2 / (2m) — E?, m is the 
effective electron mass, and Ep denotes the Fermi energy. The 
second sum in Eq. ( [Tal l is over the set W of momentum vec- 
tors with | eric | < tkJ c , ui c being the cut-off frequency. The su- 
perconducting order parameter A is assumed to have s-wave 

symmetry with A = W Skew ( c -U c +^t)- Here > w o is an 
attractive momentum-independent interaction constant. 

The superconducting system ([Tat is perturbed by a Gaus- 
sian pump pulse, which creates finite nonequilibrium quasi- 
particle distributions. In the Coulomb gauge, the optical pump 
pulse is described by a transverse vector potential 



A q (i) = A 



with full width at half maximum (FWHM) r p , amplitude Ao, 
photon frequency ui p , and photon wave vector q p . The cou- 
pling of the vector potential A q to the superconductor ( [Tal l is 



given by H em = H t 



(i) 



Hem, where 



H (2) 



2m 

er 
2m 



(lc) 



k,q,<x 

E 

k,q,cr 



'k+q,cr k,<T- 



We add to this Hamiltonian a noninteracting phonon system 
and a coupling between the phononic and electronic degrees 
of freedom. The free-phonon Hamiltonian H p h is described by 
i?ph = J2 P (bpb p + |), where b p is the annihilation op- 
erator of a phonon with wave vector p and constant frequency 
w p h. For the sake of simplicity, we restrict ourselves to a sin- 
gle branch of phonons. A generalization to several phonon 
modes is straightforward. The superconductor is coupled to 
the phononic system via the interaction 



9ph ( b -P + 6 p) C k+p, ( x < 

p,k,<7 



k , CF ' 



(Id) 



where g p ^ denotes the electron-phonon coupling constant. In 
the following we assume that the electron-phonon coupling 
strength is much smaller than the superconducting energy 
scalesr& such that the influence of the phonon subsystem on 
the superconductor becomes negligibly small. 



III. DENSITY-MATRIX FORMALISM 



Physical observables, such as the order parameter ampli- 
tude |A(t)| and the lattice displacement U(r,t), can all be 
expressed in terms of the Bogoliubov quasiparticle densities 
and the mean phonon amplitudes. Hence, we derive equa- 
tions of motion for these quantities using the framework of 
the density-matrix formalism. To this end, it is advanta- 
geous to perform a canonical Bogoliubov transformation of 
the fermionic operators, with a k = UkC k ^ + ^kcLj-j. an d 
/? k = MkC^L k ^ — WkC k .., where the coefficients and are 
time-independent and chosen such that the BCS part of the 
Hamiltonian, H sc , in the initial state, i.e., at t = ti, takes di- 
agonal form (see Appendix [A). Due to the interaction term 
iJ e . p h, Eq. ( fTdt , the equations of motion for the single-particle 
density matrices are not closed, but give rise to an infinite 
hierarchy of equations of higher-order density matrices. For 
the purpose of studying the generation of coherent phonons, it 
suffices to break this hierarchy at first order, which amounts to 
neglecting all correlations among quasiparticles and phonons. 



Thus, phonon-assisted quantities, such as (ct k a k ,6 p ), are fac- 
torized according to (a k a k ,6 p ) ~ (a k a k /)(6p). A nonvan- 
ishing (b ) corresponds to a finite displacement of the lattice 
ions. That is, the lattice displacement U(r, t) is connected to 
the coherent-phonon amplitude D p (t) = (b ) + (b_ ) via 



U(r,t) = 



2Mu) ph V 



+ip-r 



(2) 



where M is the reduced mass of the lattice ions and V the 
systems volume. 

At first order in the correlation expansion in g p h the equa- 
tion of motion for the normal quasiparticle density (a k a k ,), 
as obtained from the Heisenberg equation of motion, is given 
by 
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E ( 2k + q ) • A « h + q (4 +q « k ^) - ££,- q <4«k'- q > - ^ q K</W - ^,- q (44- q ) 



q=±q P 



-^7 E E A ^-q' ' fe q (4 +q " k ') - ^,- q (4« k '- q ) + + M i,- q (44- q ) 



2m 

q \q'=± q , 

-.9 P h E [^,-p(44'-p> + < P KA+ P ) + ik,p(4 +P «k') - ^,-p(4«k'- P ) 



(3) 



where i?k = £k(l — 2v£) + 2uk«kReA, Ck = — 2£kMk^k + 
Au k - A*^, = u k Uk+k' ± Wk«k+k', and M^ k , = 

Vkitk+k' ± w-k^k+k'- Comparing the first and the last line 
of Eq. (0), one sees that the quasiparticle-phonon interaction 
at first order in the hierarchy simply leads to a nondiagonal 
energy renormalization. The equations of motion for the re- 
maining three quasiparticle densities, (4Ac')' (44')' an< ^ 
( a kAc')' which have a similar structure, are given in Ap- 
pendixlAl 

The time dependence of the coherent-phonon amplitude 
D p (t) can be expressed in terms of a harmonic oscillator- 
type second-order differential equation (for details see Ap- 
pendixlAl 

1 D p (t) = Jp(t), (4a) 



dt 2 



with forcing term 



■Fp(*) 



2uj 



P h 



f/ph 



E Kp (< 



*k+pAc 



<44 H 



, P ((4 



k"kH 



(4 



k+p/^k/ 



(4b) 



which is purely real. Within the framework of model (Q~|), 
the equation of motion for the coherent-phonon amplitude 
Dp (£) is exact up to higher-order corrections in the correla- 
tion expansion. It is worth noting, that at the next order in 
the hierarchy (i.e., at second order in <? p h) incoherent phonons 
and quasiparticle-phonon scattering processes are generated, 
which give rise to a finite lifetime of the coherent phonons and 
which thereby lead to an exponential damping of the coherent- 
phonon oscillations. Focusing on time scales much shorter 
than the coherent phonon lifetime, we neglect in the follow- 
ing any finite lifetime effects due to quasiparticle-phonon or 
phonon-phonon scattering processes. 

Eq. (O and the corresponding equations for the other three 
quasiparticle densities (see Appendix[A} together with Eq. (HJi 
form a closed set of coupled differential equations. In SecfVl 
we solve numerically this set of equations to determine the 
temporal evolution of the order parameter amplitude |A(t)| 
and the lattice displacement U(r,t). Before doing so, we 
present in Sec. [IV] a qualitative analysis of the differential 
equation (|4|i and derive approximate solutions for different 
time scale regimes. 



IV. COHERENT PHONON GENERATION MECHANISM 

The equation of motion (O for the coherent-phonon ampli- 
tude Dp (t) resembles the equation of a forced harmonic os- 
cillator with driving force J- p (t). The forcing term -F p (i) is a 
function of the quasiparticle densities and implicitly depends 
on the optical excitation conditions, since both the normal and 
the anomalous quasiparticle densities are driven by the optical 
pump pulse. Hence, a rapid increase in the Bogoliubov quasi- 
particle distribution function due to optical excitation acts as 
a driving force for coherent-phonon oscillations. To make this 
more precise, let us express the general solution of the second 
order differential equation (HJi as 



D P (t) 



(5) 



where we assumed the following initial conditions: D p (ti) = 
and ^ D p (ti) = 0, for all p. Depending on the considered 
ordering of time scales, the time dependence of the driving 
force F p can be approximated by different functions. 

First, we focus on the regime r p <^ ta ~ r p h, where both 
the quasiparticle and phononic subsystems evolve in a nonadi- 
abatic manner and the phonon period r p h is of the same order 
of magnitude as the dynamical time scale of the order param- 
eter ta- A number of recent publications have investigated 
this regime, albeit in the absence of phonon interactions^— 
Indeed, an exact solution has been derived for the dynamics 
of a BCS superconductor after an abrupt perturbation by, e.g., 
an interaction quenchi 2 ^ 2 - 7 - i n particular, it has been shown 
that as t — > co, the absolute value of the order parameter 
|A(t)| approaches, in an oscillatory fashion, a constant value 
Aoo < |A(ii)|,i.e. 



|A(*)| 



— -p cos (2A 00 t/h- 
yt 



(6) 



where b and <fi are constants that depend on the initial state. 27 
The evolution of the normal and anomalous quasiparticle den- 
sities shows a similar oscillatory behavior with al/\/t decay. 
As it turns out, the coupling to phonons does not qualitatively 
alter this time dependence, as long as the electron-phonon in- 
teraction strength is small compared to the superconducting 
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gap amplitude. Hence, we approximate the forcing term in 
Eq. © as 

T p (t) ~ e(*)[Ap + B p cos(2A 00 t/n)/Vt] ) (7) 

with O(i) the Heaviside step-function. Inserting Eq. © 
into Eq. ©, and assuming that the phonon frequency w p h is 
close to resonance with the order parameter oscillations, i.e., 
uj& = \2Aoofh— w p h| <C w p h, we find that, to leading order in 
Wd/ojph, the coherent-phonon amplitude D p (t) is given by 



D p (t) ~ —\hr— cos{tuj vh )S 2 {tuj i ) + sin(tWph)C 2 (fcdd) 

^ph V 



(8) 



for t > 0, and where S2 and C2 denote the two Fresnel 
integrals^ In other words, the time evolution of D p (t) ex- 
hibits a beating-like phenomenon, i.e., D p (t) oscillates with 
frequency cj p h and an amplitude that is modulated by the 
Fresnel integrals (cf. Figs. Q] and ©. Exactly at resonance, 
hoj p h = 2Aoo, the coherent-phonon amplitude takes the form 



D P (t) 



w P h 



cos(wphi)] 

^ ^=sin[ Wph (i-20], (9) 
w ph Jo 2V*' 



for i > 0. As f increases, the second term quickly dominates 
in the above expression and, hence, the amplitude of the os- 
cillations in D p (t) grows like \ft. This is in excellent agree- 
ment with the numerical simulations presented in Sec.[Vl(cf. 
Fig.EJ. 

Second, we consider the regime t p ,ta <C r p h, where the 
Bogoliubov quasiparticle oscillations average out on the time 
scale of the phonons. In this case, provided that A p is not 
too small compared to B p in Eq. ©, the forcing term T p (t) 
can be approximated by J- p (t) ~ A p Q(t). Inserting this into 
Eq. © yields for the coherent-phonon amplitude D p (t) 



D P (t) 



-4 r 



[1 - COs(Wph^)] 



for t > 0. 



(10) 



Thus, the phonon oscillations are cosine-like, with the ex- 
trema lying at integer and half-integer multiples of the phonon 
period r p h. The amplitude of the oscillations increases with 
decreasing phonon frequency w p h. Again, we find good agree- 
ment with the numerical results of Sec. |V](cf. Fig. [3]). Note 
that when A p /B p becomes sufficiently small in Eq. ©, then 
there appear fast oscillations with frequency 2A oc /fi, super- 
imposed on the slow oscillations of Eq. dTob (cf. solid black 
and dotted red curves in Fig. [3j- 



V. NUMERICAL SIMULATIONS 

In this section, we numerically solve the closed set of 
equations of motion, Eqs. ©, ©, ( lAU and (1A2K both for 
the quasiparticle densities and the mean phonon amplitudes. 
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FIG. 1: (Color online) Panel (a): Numerical simulation of the lattice 
displacement (7(0, t) versus time t for three different pulse widths 
r p = 15 ps (solid black), 5 ps (dashed blue), and 0.05 ps (dotted 
red). The integrated pump pulse intensity for each trace is chosen 
such that 2AOQ = 1.7 meV. Here we take dj p h = 2.0 meV/fi. Panels 
(b) and (c) show the spectral distribution of £/(0, t) and the temporal 
evolution of | A(t)j, respectively, for the same parameter values as in 
panel (a). In panel (c), the curves are vertically shifted by multiples 
of 0.1 meV. 



From these quantities, the temporal evolution of the lattice 
displacement U(r,t) and of the order parameter amplitude 
|A(f)| is readily computed. Inspection of Eq. © shows that 
the largest entries in the quasiparticle density matrices are 
those confined to a band centered around the diagonal. That 
is, off-diagonal entries, such as e.g. (a^k+nqp) > are °f order 
|Ao| n . Hence, for sufficiently small |An|, the off-diagonal 
elements decrease rapidly as n increases. To reduce the com- 
putational effort we therefore set all off-diagonal entries with 
n > 4 to zero and, furthermore, restrict ourselves to a one- 
dimensional wire geometry. It is important to note, however, 
that the phenomena discussed in this paper are qualitatively 
independent on the dimensionality of the system, as can be 
seen from the analytical analysis given in Sec. [IV] In fact, for a 
related model it has been shown that (quasi-)one-dimensional 
simulations provide a good approximation for two- and three- 
dimensional superconductors. 2 - 

For the numerical computations we use the following ma- 
terial parameters*^ superconducting gap in the initial state 
A(tj) = 1.35 meV, cut-off energy hu c = 8.3 meV, Fermi 
energy Ep = 9479 meV, effective electron mass m = 1.9mo, 
with tuq the free electron mass. The optical pump pulse is 
centered at t = and has a central energy of hu) p = 3 meV, 
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FIG. 2: (Color online) Panel (a): Lattice displacement U (0, t) as a 
function of time t for resonant (dotted red: Aqt p = 415, 2Aoo = 
2.0 meV) and off-resonant conditions (solid black: Aqt v = 254, 



2A C 



2.3 meV; dashed blue: A&t v = 561, 2A C 



1.7 meV). 



Here we take r p = 0.5 ps and cjph = 2.0 meV/ft. The gray curve 
shows the \/t dependence predicted by Eq. Panel (b): Spectral 
distribution of the coherent phonon oscillations for the same param- 
eters as in panel (a). Panel (c): Asymptotic value of the order param- 
eter Aoo versus integrated pump pulse intensity for seven different 
pulse widths Tp, 



which is of the same order, but slightly larger than 2A(ij). 
As the initial state for the simulations we choose the equi- 
librium BCS ground state at zero temperature. We find that 
the lattice displacement U(r, t) has a quite weak dependence 
on position, exhibiting oscillations in time at all values of r. 
The frequency of these oscillations does not depend on posi- 
tion, but their amplitude varies weakly with r. Therefore, we 
choose to present in all the figures only the displacement for 
r = 0. A detailed study of the r dependence of U(r,t) will be 
presented elsewhere. 

In the following, we adjust the pump pulse length r p , the 
phonon frequency w p h, as well as the integrated pump pulse 
intensity AqT p , to explore different time scale regimes. 

a. ta ~ T p h! We start with the most interesting case, 
namely the situation where the Bogoliubov quasipraticle os- 



cillations are close to resonance with the phonon frequency, 
see Figs. Q] and [2] In Figs. 02a) andQJc), we plot U(0,t) and 
|A(t)|, respectively, for three different pump pulse lengths 
t p = 0.05, 5, and 15 ps. The integrated pump pulse inten- 
sity for each curve in Fig.[T]is adjusted such that 2Aoo = 1.7 
meV. This ensures that the order parameter oscillations are al- 
ways close to resonance with the phonon frequency. When the 
quasiparticle subsystem is perturbed nonadiabatically (r p < 
ta, i.e. t p = 0.05 ps in Fig.Q]and r p = 0.5 ps in Fig.0, the 
quasiparticle densities build up in a coherent manner while the 
system is out of equilibrium. This leads to rapid oscillations in 
the quasiparticle densities, and hence also in the order param- 
eter, with an amplitude decaying as \j\fi and a frequency that 
is determined by the asymptotic gap value [see Fig.[TJc) 
and Eq. ©], Consequently, the coherent phonons are driven 
by a sinusoidal forcing term. Thus, whenever the order pa- 
rameter oscillations are close to resonance with the phonon 
mode (i.e., \2AoofK — cj p h| <C w p h), we observe a pronounced 
beating phenomenon [cf. Eq. {§)]. 

Most importantly, we find that the frequency of the order 
parameter oscillations can be tuned exactly to resonance by 
adjusting the integrated pump pulse intensity [cf. Fig. |2|c)]. 
This is demonstrated in Fig. [2ja), which is the main result of 
our paper. In this figure we plot the lattice displacement for 
fiw p h = 2 meV and a fixed pulse duration r p = 0.5 ps, but dif- 
ferent integrated pulse intensities. At resonance, AqT p = 415 
(2Aoo = 2 meV), the amplitude of the phonon oscillations 
shows a square-root increase with t, which is in agreement 
with Eq. (|9]l. In Figs. Q2b) and|2|b) we present the spectral 
distributions of the coherent phonon oscillations as obtained 
from the Fourier transforms of {7(0, t). The discussed behav- 
ior of {7(0, t) reflects itself in the Fourier transforms: for res- 
onant condition we observe a strong single peak at 2 meV, 
while for off-resonant condition there are two peaks, one at 
the phonon energy ficj p h = 2 meV and the other at the fre- 
quency of the order parameter oscillations, i.e. at 2Aoo = 1.7 
and 2.3 meV, respectively. 

In the regime, where the quasiparticles evolve adiabati- 
cally (t p > ta, i.e. t p = 5 ps and 15 ps in Fig. Q]), the 
pump pulse drives only the normal quasiparticle densities, 



and 



and (/3tj3 k ,), but leaves the anomalous ones, (a k /3 k , 



t k /? k ,), mostly unaffected^ Thus, the instantaneous 



value of the gap is almost fully determined at all times by 
the quasiparticle occupations, (a k a k ,) and (/3 k /? k ,), and the 
gap amplitude decreases monotonically from its initial equi- 
librium value A(tj) to its final value A M [black solid curve 
in Fig. 0I C )]- As it turns out, in this situation the coherent 
phonons are still driven by a sinusoidal forcing term of the 
form ((TJl, albeit with a much smaller amplitude. As a conse- 
quence, U (r, t) still exhibits a beating phenomenon, but has a 
considerably smaller magnitude than in the nonadiabatic case. 
Deep inside the adiabatic regime (t p 3> ta) the coherent 
phonon oscillations eventually vanish completely. 

In passing, let us also comment on the dependence of 
on the integrated pump pulse intensity AqT p , which is shown 
in Fig. |2tc) for seven different pulse widths t p . The asymp- 
totic gap value A^ is linear at small AqT p for all t p , but de- 
viates from this linear behavior at higher integrated intensi- 
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FIG. 3: (Color online) Lattice displacement U (0, t) as a function of 
time t for three different pump pulse lengths r p = 0.05, 2, and 10 ps 
and two different phonon frequencies u p h = 0.1 and 0.05 meV//i, 
corresponding to r p h = 41 and 83 ps, respectively. The integrated 
pump pulse intensity for each curve is adjusted such that 2Aoo = 1.7 
meV. The curves with r p = 2 and 10 ps have been multiplied by the 
constant factor 50 and 400, respecitvely. 
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FIG. 4: (Color online) Lattice displacement U (0, t) versus time t for 
three different integrated pump pulse intensities Aqt p = 129, 343, 
and 969. These integrated intensities correspond to 2AOO = 1.7, 
2.0, and 2.3, respectively. Here we choose r p = 15 ps and uj p h = 4 
meV/fi. The gray trace depicts the Gaussian time dependence of 
the pulse envelope. The inset shows the spectral distribution of the 
coherent phonon oscillations for the same parameters as in the main 
panel. 



ties. While the curves corresponding to short pump pulses 
(r p = 0.05 and 0.5 ps) exhibit a downward bend, those with 
longer pulse widths (r„ = 2.0, 5.0 and 20.0 ps) flatten with 
increasing AqT p . (The curves with r p = 1.0 and 1.38 ps 
lie in between these two regimes, showing first a downward 
and then an upward bend.) The downward bend is due to a 
quadratic term in the AqT p dependence resulting from two- 
photon processes. The flattening, on the other hand, occurs 
because long pump pulses, with r p ^> 2tt/uj p ~ 1.38 ps, 
create sharp and narrow peaks in the quasiparticle distribu- 
tions, which, for sufficiently high intensities, leads to satura- 
tion due to Pauli blocking. 28 - 1 We observe that the integrated 
intensity above which Pauli blocking sets in decreases with 
increasing r p . 

b. ta, t p <C t p i,: In Fig. [3]we show the lattice displace- 
ments U(0,t) induced by optical pump pulses with pulse 
lengths T p = 0.05, 2 and 10 ps for two different phonon 
energies fiw p h = 0.1, and 0.05 meV. This parameter choice 
corresponds to the case, where the rapid oscillations in the 
quasiparticle subsystem average out on the time scale of the 
phonons. Since t p <C T p h, the phonons are perturbed by 
an almost instantaneous change in quasiparticle occupations, 
which leads to cosine-like coherent phonon oscillations with 
frequency w p h and an amplitude that increases with decreas- 
ing phonon frequency w p h [cf. Eq. ([TOiil. As we go from 
the regime where the Bogoliubov quasiparticles are perturbed 
nonadiabaticallly (t p < ta, i.e. t p = 0.05 ps in Fig. |3j to 
the regime where the quasiparticles are perturbed adiabati- 
cally (t p > ta, i.e. t p = 10 ps in Fig. |3]l the amplitude 
of the oscillations decreases quickly. For sufficiently long 
pump pulses, eventually there appear fast oscillations with fre- 
quency 2Aoo/?i on top of the slow oscillations with frequency 
w p h (t p = 10 ps in Fig.E). 



c. T p h < ta <C t p : Lastly, we consider the case where 
both the quasiparticles and the phonons are perturbed in an 
(almost) adiabatic fashion and the coherent phonon oscilla- 
tions are off resonance. In Fig. @] the time dependence of 
U(0, t) is shown for a pump pulse with length t p = 15 ps and 
three different integrated pump pulse intensities. The phonon 
energy is chosen to be ftw p h = 4 meV, which is larger than 
2|A(fj)| and hence far away from resonance. In spite of 
the almost adiabatic evolution of the system on the phonon 
time scale, coherent phonons are still being generated, albeit 
with a much smaller amplitude than in Figs.[T][3] Remarkably, 
the Fourier spectrum of the coherent phonon oscillations (in- 
set in Fig. HI does not only show contributions at 2Aoo and 
hoj p h = 4 meV, but also a third peak at 3 meV, which is iden- 
tical to the pump photon energy. The latter contribution is 
caused by large transient oscillations occurring in the time 
interval ~ [—Tp, +t p ] during which the pump pulse acts on 
the system. If the pump pulse frequency w p is chosen close 
to w p h, the coherent phonon oscillations show some enhance- 
ment, i.e., the coherent-phonon amplitude increases (almost) 
monotonically until t ~ ps, and then remains constant at its 
peak value even after the pump pulse has passed (not shown). 



VI. CONCLUSIONS 

In this paper we have presented a theoretical investigation 
of the generation of coherent phonons in a superconductor by 
ultrafast laser pulses. Using the density-matrix formalism, we 
have performed numerical simulations of the nonequilibrium 
dynamics of a BCS s-wave superconductor coupled to a sin- 
gle branch of optical phonons. Based on both numerical and 
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analytical arguments, we have shown that sudden changes in 
the Bogoliubov quasiparticle densities created by the optical 
pump pulse act as a driving force for coherent phonon oscil- 
lations. For ultrafast laser excitations, the superconductor ex- 
hibits a nonadiabatic coherent dynamics that is characterized 
by rapid order parameter oscillations. We have found that the 
creation of coherent phonons is resonantly enhanced when the 
period of these gap oscillations coincides with the phonon pe- 
riod. In a pump probe experiment this resonance condition 
can be achieved by tuning the frequency of the gap oscilla- 
tions via a change in the integrated pump pulse intensity (see 
Fig.©. 

The resonant coherent phonon generation mechanism dis- 
cussed in this paper applies in principle to any BCS-type su- 
perconductor that has an optical phonon with phonon energy 
of the same order as the superconducting gap. One interest- 
ing class of examples are superconductors which are close to 
a structural transition that is driven by a soft optical phonon, 
i.e., e.g., CaC u 40 ' 41 or CaAlSi. 42,4? The coherent phonon os- 
cillations are experimentally observable, for example, as peri- 
odic modulations in time-resolved reflectivity measurements. 



Driving a superconductor into the regime of nonadiabatic co- 
herent dynamics requires ultrashort laser pulses with frequen- 
cies of the order of the superconducting gap, i.e., in the ter- 
ahertz regime. With the recent advent of ultrafast terahertz 
sources, 44 we hope that it will be soon possible to perform 
time-resolved measurements on superconductors in the nona- 
diabatic regime and to test our theoretical predictions. The 
experimental observation of the discussed resonant coherent 
phonon oscillations would not only be interesting in itself, it 
could potentially also give useful information about the gap 
symmetry and the pairing mechanism of the superconductor. 
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Appendix A: Equations of motion 



In this Appendix, we give the equations of motion for the quasiparticle density (44') and (4 Ac')' and the mean phonon 
amplitude (&_). In deriving these differential equations, we use the fact that both u k and u k are real and time-independent, with 

As explained in the main text, we 



u k = ^1/2(1 + e k /£ k ) and v k = ^1/2(1 - e k /£ k ), and where E k = ^ + | A(t,)| 2 . 
neglect terms of second- or higher-order in the correlation expansions in g p h and decouple phonon-assisted quantities according 

to, e.g., (44'^p) — (44') (^p)- By use of Heisenberg's equation of motion, we find that the time dependence of (44') is 
described by the following differential equation 



d_ 

"dt 
eh 

2m 



0R k + i?kO (44') + c u> (4«ic> + c t (4 Ac) - s- 



3k', k 



J2 (2k + q) • A q U+ q (4 +q 4> - 4,- q (44- q > + ^,- q (4« k '- q ) + M k - q (/3 k+q 4) 



q=±q P 



2m 



E E a 



q - q ' • A q / 



^k, q (4 +q 4> + £ k ',- q (44- q > + m+ _ q <4«*_ q ) - M+ q (/3 k+q/ 3 k ,) 

q \ q '=± qt 

-.9phE^K4 (<4a+p> -<w<) +£ k ',- P (44-p> + ^k,p(4 +P 4'> + m+ _ p <4« k ,_ p 
p 



(Al) 



where i? k , C k , -L k , and M k are defined in Sec. HIT] Note that the equation of motion for (a k /? k ,) can be obtained straightfor- 
wardly from Eq. (lAlb by complex conjugation. The equation of motion for (4 Ac') rea ds 



* h j t (4 Ac) = {R* - R k ) (4 A") + c k '(44) + c k K/3 k ,) 



eh 
2m 



£ (2k - q) • A q fe,_ q {4- A> - 4, q (4Ac +q } - M k -_ q (a k - q Av) - ^, q (4 +q 4) 



q=±q P 



2m 



E E a 

q \ q '=± qp 



q - q ' • A q / 



£k,- q (4- q Ac> - ^ k ', q (4/3k' +q > + M+ q <4, +q 4) + A/+„ q (a k _ q /3 k /) 



E [K-^-M + A4+ P (4 + p4> - k>,Mpk+p) + L k,- P (4-p/v; 



(A2) 
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We also present here the equation of motion for the mean phonon ampliuted (b ), which is derived in a similar fashion as 
Eq. dAD, 

ihj t {b p ) = fi^ ph (6 p ) - (A3) 

where T v (t) is defined in Eq. ( f4bb - Again, we note that the equation of motion for (6_ p ) is related to the one for (6 ) by complex 

conjugation. Adding the equations for (bJ) and (b_ p ) and taking a time derivative, one derives the equation of motion for the 
coherent-phonon amplitude D p (t), Eq. ©. 

As can be seen from Eqs. ||3}, (lAU . and ( IA2I ). any element in the quasiparticle density matrices, such as e.g. (aj c o: k+Ilq ) 
is strongly coupled only to those elements with indices in the subspace (k + Zq p , k + lq p + mq p ), where I, m, and n are 
integers. Elements with indices in different subspaces are only weakly coupled through the superconducting order parameter. 
For the one-dimensional simulations we therefore discretized the momentum space by a one-dimensional grid with mesh size 
|q p |. Furthermore, we approximate (2k + q) ■ A q in the second line of Eqs. ||3}, ( lAU . and ( IA21 > by 2kp • A q , which is justified 
since the photon wave vector is much smaller than the Fermi momentum |kp|. We use a standard fourth-order Runge-Kutta 
technique to integrate the equations of motions given by Eqs. (|3), ©, dAU and (IA2t . 
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